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The Stochastic Variational Method (SVM) is used to show that the effective mass model correctly 
estimates the binding energies of excitons and trions, but fails to predict the experimental binding 
energy of the biexciton. Using high-accuracy variational calculations, it is demonstrated that the 
biexciton binding energy in transition metal dichalcogenides is smaller than the trion binding energy, 
contradicting experimental findings. It is also shown that an excited state of the biexciton is in very 
good agreement with experimental data. This excited state corresponds to an hole attached to a 
negative trion and may be a possible resolution of the discrepancy between theory and experiment. 


Since the isolation of graphene, much attention has 
been directed toward two-dimensional (2D) materials and 
their extraordinary potential within nanoelectronics and 
photoelectronics [1-9]. Of specific interest in recent re¬ 
search has been the development of a graphene-like 2D 
material featuring a sizable band gap, allowing it to serve 
as a 2D semiconductor. In particular, transition metal 
dichalcogenides (TMDs) have been investigated as pop¬ 
ular candidates. An observed consequence of reduced 
dimensionality and weak dielectric screening in such ma¬ 
terials is a strong electrostatic interaction allowing the 
formation of bound electron-hole pairs (excitons) with 
binding energies on the scale of several hundred meV in 
TMDs such as WSe 2 [9]. This behavior is not observed 
within 3D bulk counterparts and is thus notably unique 
to few-layer materials. Charged excitons (trions) have 
also been observed with surprisingly large binding ener¬ 
gies in M 0 S 2 (20 meV) [10], WS 2 (45 meV) [11], and 
WSe 2 (30 meV) [12]. Such large exciton and trion bind¬ 
ing energies suggest the existence of bound exciton pairs 
(biexcitons) in monolayer TMDs, which have indeed been 
experimentally found in WSe 2 (52 meV) [13] and WS 2 
(65 meV) [14]. 

There is a substantial need to understand these exci- 
tonic structures in 2D materials in order to character¬ 
ize their electrical and optical responses and fully assess 
their potential functionality. Effective-mass models are 
often used to model electron-hole systems [15-17] and 
have been widely applied to excitonic systems in 2D ma¬ 
terials [18-22]. These models have been successful in the 
calculation of binding energies and geometries (extended 
wavefunctions) of excitons and trions in TMDs. 

In this work, we show that while the effective mass 
model correctly estimates the binding energies of exci¬ 
tons and trions, it fails to predict the binding energy of 
the biexciton. Using high-accuracy variational calcula¬ 
tions, we demonstrate that the biexciton binding energy 
predicted in TMDs by the effective mass model is smaller 
than the corresponding trion binding energy, contradict¬ 
ing the aforementioned experimental findings. However, 
we also show that the binding energy of an excited state 
of the biexciton closely agrees with experimental data. 


This state corresponds to an electron attached to a posi¬ 
tive trion and may resolve the discrepancy between the¬ 
ory and experiment. 

The computational technique employed in this paper 
is the Stochastic Variational Method applied to exciton 
(X = eh ), trion (X - = eeh), and biexciton (X 2 = eehh) 
systems using the explicitly correlated Gaussian (ECG) 
basis. This basis is known to produce high-accuracy 
binding energies (8-10 decimal digits) when applied to 
similar few-particle systems, including H 2 [23], Hj, and 
the positronium molecule PS 2 [24, 25]. 

The nonlinear parameters of these explicitly correlated 
Gaussians are optimized using stochastic variation, a pro¬ 
cedure which uses random trial and error to iteratively 
improve the quality of the ECG-basis representation of 
the desired wavefunction. This method is known to be 
well-suited to the description of both ground and excited 
states of excitonic structures with up to five particles, 
such as positively charged biexcitons Xj (eehhh) [16, 26]. 
We refer the reader to Ref. [27] for a through review of 
the applications of the ECG basis in various problems. 

The nonrelativistic Hamiltonian of an excitonic N- 
particle system is given by 



where iq, m^, and qi are respectively the 2D position 
vector, effective mass, and charge of the ith particle, 
rij = \ri — Tj\, and r s is the screening length of the 
medium. The 2D screened electrostatic interaction po¬ 
tential is then given by 

V(r) = ^[H 0 (r)-Y 0 (r)} (2) 

where Hq and lo are the Struve function and Bessel 
function of the second kind, respectively. This poten¬ 
tial differs substantially from the usual 1/r Coulomb po¬ 
tential, exhibiting nonlocal macroscopic screening which 
arises in 2D systems [28]. For small distances (r —» 0) 
the potential exhibits logarithmic divergence, while for 
large distances (r —>> oo) it asymptotically falls off as a 
1/r Coulomb potential. The 2D layer polarizability \ 2 D 
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determines the length scale r s = 2 ttx 2 D that separates 
these two regimes. 

As a variational trial function, we adopt a 2D form of 
the correlated Gaussians [27, 29] 

^ X A aPi' | 

( 3 ) 

where 

tm(p) = (Px + iPy)" 1 ■ (4) 

Here, pi denotes the ith relative (Jacobi) coordinate of 
the system, nrii denotes the magnetic quantum number 
of the ith particle, and M — m i + m 2 + • • • + mjv- This 
function is coupled with the spin function Xsm s t° fc> rm 
the trial function. 

The material-dependent parameters (effective masses 
and screening lengths) adopted in this paper are those 
given in Ref. [20]. These parameters are based on ab 
initio calculations and produce exciton binding energies 
in close agreement with Bethe-Salpeter calculations. For 
simplicity, we fix the electron-hole mass ratio a = m e /mh 
to unity, and consider only negative trions (eeh) in our 
calculations. (No loss of generality is incurred, as when 
cr = l, charge conjugation symmetry guarantees that the 
energies produced are equivalent to those of positive tri¬ 
ons.) Indeed, experimental observations find that the 
binding energies of positive and negative trions are nearly 
equal [5], validating our choice of a = 1 as a good ap¬ 
proximation. Even so, we point out that slight variations 
in a do not affect our qualitative results. 

The calculated and experimental binding energies of 
excitons are reported in Table I. Note that the ex¬ 
perimental binding energies are in the 500 meV range 
[3, 6, 8, 9] and the parameters given in Ref. [20] give good 
overall agreement with experiment. More accurate agree¬ 
ment is not pursued, partly due to the presence of uncer¬ 
tainties in the experiments (e.g.substrate dependence), 
and partly because the binding energies of the trions and 
biexcitons are not very sensitive to the exciton binding 
energy. The present SVM approach reproduces the en¬ 
ergies of Ref. [20] for excitons and improves the results 
of Ref. [20] by a few meV for trions. This is reason¬ 
able because we use a variational ansatz with 100 basis 
functions, while in Ref. [20] only a single trial function is 
used. 

The experimental observation of biexcitons in mono- 
layer WSe 2 reported in Ref. [13] is accompanied by a 
variational calculation estimating the binding energy to 
be 37 meV with a one-term variational trial function. 
This is lower than the measured value of 52 meV reported 
in the same paper. The authors expect that a more ac¬ 
curate approach would reconcile theory and experiment, 
but our calculation shows that this is not the case. A 
larger basis provides improved estimates of both exciton 


4>m( r) = A 



and biexciton ground-state energies, giving a biexciton 
binding energy even smaller than the single term predic¬ 
tion in Ref. [13] (see Table I). 

The underbinding of the biexciton is not unique to 
WSe 2 - We consistently observe that for typical screening 
lengths occurring in M 0 S 2 , MoSe 2 , WS 2 , and WSe 2 , the 
biexciton binding energy is always significantly smaller 
than the trion binding energy (see Table I). This stands 
in contrast to experiment, where trion binding energies 
are usually found in the 20-30 meV range, and biexciton 
binding energies are found between 50 and 70 meV. A 
general trend appearing in our calculations is that above 
a certain screening length (r s > 5 a.u.) the biexciton 
is less strongly bound than the trion, while for smaller 
screening lengths (r s < 5 a.u.) the biexciton becomes 
more strongly bound. However, even this adjustment in 
screening length cannot be made to agree with experi¬ 
ment, since in this range the exciton binding energy dif¬ 
fers significantly from the observed value. (The exciton 
binding energy increases to 2.5 eV, while the trion bind¬ 
ing energy becomes 250 meV). Varying the electron-hole 
mass ratio a would slightly change the results but again 
would not help to resolve the disagreement. 

In addition to its ground state, the biexciton also has 
three bound excited states [16, 30, 31]. Two excited 
states exist with L = 0 and positive parity, and one exists 
with L = 1 and negative parity [30]. One of the L = 0 
excited states is bound due to charge inversion symme¬ 
try, but this is only valid in the case of equal electron and 
hole masses. Thus, we have not considered this state in 
our calculations. The L = 0 excited state is bound with 
respect to the X(1S)+ X(2 S) threshold (where X(nL) de¬ 
notes the nth exciton state with orbital momentum L), 
and the L = 1 state is bound with respect to the X(15) + 
X(2 P) channel. (Note that while the usual 1/r Coulomb 
potential gives equal energies for the 2 S and 2 P exciton 
states, this degeneracy is broken by the screened 2D po¬ 
tential under present consideration.) These excited states 
cannot autodissociate to the X(15)+X(15) threshold, as 
symmetry considerations [30-32] force this channel to be 
closed. 

Our calculations show that the energy of the L = 1 ex¬ 
cited state is lower than that of the L = 0 excited state, 
while simultaneously the threshold of the L = 1 state 
is higher than the threshold for the L = 0 state. This 
results in binding energies of order of 250-450 meV (see 
Table I), much larger than experimental results. More¬ 
over, this L = 1 state is much less likely to be formed 
than the L = 0 state. 

To investigate the structure of the excited state, we 
studied the pair correlation function 
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where p and q stand for electrons or holes, and the sum is 
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TABLE I. Comparison of experimental and theoretical TMD excitonic binding energies (meV). 



System 

M0S 2 

MoSe2 

WS 2 

WSe 2 

Experiment 

X 

500 [33], 570 [34] 

550[35] 

320 [36], 700 [3] 

370[9] 

Theory 

X 

555 

480 

523 

470 

Experiment 

X“ 

18±1.5[10] 

30 [5, 37] 

30 [14], 45[11] 

30 [8, 12] 

Theory 

X“ 

34 

28 

34 

30 

Experiment 

x 2 

70 [33] 

- 

65[14] 

52 [13] 

Theory 

x 2 

22 

18 

24 

20 

Theory 

x^ (L = 0) 

69 

58 

67 

59 

Theory 

X^ (L = 1) 

460 

430 

360 

240 



P 



FIG. 1. (Color online) Electron-electron (a) hole-hole (b), 
electron-hole (c) correlation functions for the WSe 2 trion 
(solid line) ground-state biexciton (dashed line), and L — 0 
excited biexciton (dotted line). Atomic units used. 


taken only over corresponding pairs. Fig. 1(a) shows the 
electron-electron (ee) correlation functions. Those of X - 
and X 2 are quite different. By adding a hole to X - and 
thus forming a singlet state with the hole of the trion, 
the resulting X 2 becomes more compact than the trion 



FIG. 2. (Color online) Schematic view of the excited biexciton 
as a system of a trion core (red) and a long tail representing 
the hole (blue). 


due to the extra attraction. The ee correlation function 
of the X - and X^L = 0) are, however, almost identical, 
indicating that the hole (which is in this case coupled as a 
triplet with the hole of the trion) is situated somewhere 
far away from the trion and thus does not disturb the 
structure of the trion. 

Fig. 1(b) shows the hole-hole (hh) correlation func¬ 
tions. These are different for the ground and excited state 
of the biexciton, as expected from the previous discus¬ 
sion. The electron-hole (eh) correlation functions, shown 
in Fig. 1(c), are very similar in X - and X 2 ; the elec¬ 
tron and hole take up the most energetically favorable 
positions in the trion, and because there is no symmetry 
restriction, the same occurs in the ground state biexci¬ 
ton. In X^(L = 0), however, the spin triplet nature of 
the hole wavefunction restricts the available spatial posi¬ 
tions. The peak of the X^L = 0) eh correlation function 
is at the same position but roughly half the amplitude of 
that of X - and X 2 , indicating that one of the holes exists 
in a trion-like structure. The tail of the X^L = 0) corre¬ 
lation function extends far beyond the other correlation 
functions, showing that a hole exists somewhere outside 
of X - and corroborating the X - + h structure shown in 
Fig. 2. These extended correlation functions also sug¬ 
gest that simple one-term wavefunctions are unlikely to 
provide accurate descriptions of these structures. 

In summary, by using high-accuracy variational cal- 
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dilations, we have shown that in models aiding exper¬ 
iments, the binding energy of biexcitons in transition 
metal dichalcogenides is less than half of the experimental 
value. An excited biexciton, however, is a possible can¬ 
didate for the experimentally observed state. Investiga¬ 
tions of other effects, e.g. surface-bound or defect-bound 
excitons or excitonic complexes are necessary before the 
experimental data can be unambiguously assigned to a 
new biexciton formation. 
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